Exploiting Quantum Parallelism 
To Simulate Quantum Random Many-Body Systems 
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We present an algorithm that exploits quantum parallelism to simulate randomness in a quantum 
system. In our scheme, all possible realizations of the random parameters are encoded quantum 
mechanically in a superposition state of an auxiliary system. We show how our algorithm allows for 
the efficient simulation of dynamics of quantum random spin chains with known numerical methods. 
We also propose an experimental realization based on atoms in optical lattices in which disorder 
could be simulated in parallel and in a controlled way through the interaction with another atomic 
species. 
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One of the most remarkable features of quantum me- 
chanics is that it allows for the creation of superposition 
states, making the ambitious dream of performing many 
tasks at the same time a real possibility. This extraordi- 
nary concession lies at the heart of quantum computation 
and is the basis of all quantum algorithms developed so 
far pj. In this letter we present a new algorithm that 
exploits this quantum parallelism to simulate in parallel 
many different evolutions of a quantum system. Our mo- 
tivation is the simulation of physical systems whose un- 
derstanding requires to study their behavior under many 
different Hamiltonians. An important example of such 
systems are quantum random systems (QRS) 2]. For 
them certain parameters of the Hamiltonian (e.g., inter- 
action couplings, potential strengths) are random (clas- 
sical) variables. Therefore the exact simulation of their 
dynamics requires to perform many evolutions, one for 
each realization of the set of random variables. QRS have 
captured a lot of attention in the last decades 0, 0, • 
The presence of randomness can dramatically change the 
behavior of quantum many-body systems, leading to fas- 
cinating phenomena ;2]. Moreover, the answer to puz- 
zles as the unusual transport properties of high tempera- 
ture superconductor materials is inextricably tied to the 
understanding of phase transitions and transport in the 
presence of disorder p| . On the theoretical side, the un- 
derstanding of QRS is hindered by the fact that the num- 
ber of required simulations for an exact calculation scales 
exponentially with the number of random parameters Q . 
On the experimental side, one of the big challenges is the 
creation of randomness in a controlled way. Here, atomic 
systems in optical lattices 0, Q i highly versatile and con- 
trollable, are one of the most promising candidates. 

In this letter we present an algorithm that allows to 
simulate dynamics (and ground state properties) of a 
QRS within one single time evolution in which the system 
is put into interaction with an auxiliary system. The key 
idea is that all possible realizations of the set of random 
classical parameters are encoded quantum mechanically 
in a superposition state of the auxiliary system. Choosing 



the interaction with the ancilla in the appropriate way, 
all possible quantum evolutions of the QRS are simulated 
in parallel. As a particular case, adiabatic evolution with 
the ancilla can simulate at once all possible ground states 
of the QRS. As one of the main results of this work, our 
algorithm establishes an exact mapping between a QRS 
and a certain interacting (non-random) system. This 
equivalence opens a new path both to the numerical and 
experimental simulation of QRS. On the numerical side, 
it allows for the efficient simulation of QRS within the 
framework of numerical methods that simulate the cor- 
responding interacting systems efficiently. For example, 
for the case of a quantum random spin chain we will show 
that the problem is mapped onto the simulation of the 
time evolution of a one-dimensional (ID) lattice system. 
Recently, numerical methods inspired in density matrix 
renormalization group (DMRG) |9| and matrix-product- 
state (MPS) techniques have been developed to effi- 
ciently simulate the time evolution of ID lattice systems 
[llL Ha . Il3| . Here we will show how to implement the ef- 
ficient simulation of quantum random spin chains within 
the methods introduced in [T^ . On the experimental side 
our scheme opens the possibility of simulating random- 
ness in parallel through the interaction with an auxiliary 
quantum system. Moreover, conversely, it allows to in- 
terpret current interacting experimental schemes as po- 
tential simulators of certain random equivalent problems. 
We propose an experimental scheme in which a variety of 
disordered phases as quantum glasses 0, 0] or Ander- 
son insulator phases [Tg could be simulated in current 
experiments with optical lattices \PA ITs| . 

The algorithm. Let us consider a quantum system with 
Hilbert space H that evolves accordingly to a Hamil- 
tonian H(r±, . . . ,r n ) where n, . . . ,r n are random vari- 
ables that take values within a finite discrete set, re € 
Tg = {\\ , . . . A~j }, with a probability distribution given 
by p(r\, . . . , r n ). In order to simulate exactly the dy- 
namics of such a system one would need to perform 
Y[e—i m i simulations, one per each possible realization 
of the set of random variables r = (ri,...r n ). For 
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each realization the system evolves to a different state 
\ipr{t)) = e~ lH(rS>t / h \ipo}, where \tp ) is the initial state. 
Given this set of evolved states and a physical observable 
O, one is typically interested in the average of the expec- 
tation values of that observable in the different evolved 
states, that is, in quantities of the form: 

((6(t))) :=J2p(r)(Mt)\0\Mt))- (1) 

r 

On the following we describe an algorithm that allows 
to simulate in parallel all possible time evolutions of the 
random system described above. We consider an aux- 
iliary system with Hilbert space H a and a Hamiltonian 
acting onH<3H a of the form H = H(R\, . . . , R n ), where 
Rl, . . . , R n are operators that act in Ti. a , commute with 
each other and have spectra Ti, . . . , r„. Note that we 
have replaced the set of random variables r by a set of 
quantum operators R with the same spectra. The algo- 
rithm works as follows. 1) Initialization. Let us prepare 
the auxiliary system in an initial superposition state of 
the form: 

r 

where the states |r) are simultaneous eigenstates of the 
set of operators R, with Re\r) = re\r). Each state |r) 
is therefore in one to one correspondence with one real- 
ization of the the set of random variables r, its weight 
in the superposition state @ being equal to the prob- 
ability with which the corresponding realization occurs 
for the random system. 2) Evolution. We evolve the ini- 
tial state of the composite system \ip ) ® \ip a ) under the 
Hamiltonian H . The evolved state is 

\*{t)) = Y,VW)mt))®\r)- (3) 

r 

This superposition state contains the complete set of 
evolved states we are interested in. 3) Read-out. In order 
to obtain the quantities Q we just need to measure the 
observable O <g> 1, 

<¥(t)|d®l|*(t)> = <<d(t)». (4) 

The algorithm above allows us, in particular, to obtain 
the averaged properties of a random system over the 
collection of all possible ground states. Let us assume 
that the interaction between the system and the ancilla 
is introduced adiabatically, so that the Hamiltonian is 
now H(t) — H(/3(t)R), where (3(t) is a slowly varying 
function of time with /3(0) = 0, /3(T) = 1, T being the 
time duration of the evolution. If the system is prepared 
in the ground state of the Hamiltonian H(0), the algo- 
rithm above will simulate in parallel all possible adiabatic 
paths, so that the composite superposition state (|3J| will 
contain all possible ground states of the random system 

E3- 



Additionally, the scheme above can be easily extended 
for the computation of other moments of the distribu- 
tion of physical observables (higher than Q), which are 
sometimes important in the understanding of QRS Q. 
For example, quantities like {(O 2 ) — (O) 2 ), can be com- 
puted by using an additional copy of the system [2(J • 

Numerical implementation. The algorithm described 
above reduces the simulation of a quantum random sys- 
tem to the simulation of an equivalent non-random inter- 
acting problem. This exact mapping allows us to inte- 
grate the simulation of randomness in quantum systems 
within the framework of numerical methods that are able 
to efficiently simulate the corresponding interacting prob- 
lem. As an illustrative example we consider the case of 
a ID spin s = 1/2 system with random local magnetic 
field. The Hamiltonian of the system is: 

N 

H(b 1 ,...,b N ) = H +Bj2hSt, (5) 

where Hq is a short range interaction Hamiltonian, b = 
(pi, . . . , 6/v) is aset of classical random variables that take 
values {1/2,-1/2} with probability distribution p(h). 
Following the algorithm above the 2 N simulations re- 
quired for the exact simulation of the dynamics (or the 
ground-state properties) of this random problem can be 
simulated in parallel as follows. We consider an aux- 
iliary ID spin cr = l/2 system. We prepare this an- 
cilla in the initial state \i(j a ) = X^b a b|b), where the 
states |b) have all z components of the N spins well de- 
fined, af |b) = 6f|b), and «b = \/p(b). The entangled 
properties of the state of the ancilla reflect the classi- 
cal correlations among the random variables. For ex- 
ample, for a uniform distribution of the random field, 
p(b) = 1/2 , the state of the ancilla is just a product 
state, \ip a ) oc (| f) + | 1))® N ■ We evolve the system and 
the ancilla under the interaction Hamiltonian 

H = H(a{, ...,di)=H + pJ2 d ^t- ( 6 ) 

i 

Here, (3 — B if we want to simulate dynamics under 
Hamiltonian J5J, and (3 is a slowly varying function of 
time with (3(0) = and (3(T) = B for the simulation of 
the ground state properties. We have then reduced the 
simulation of the random problem to that of the time evo- 
lution of two coupled spin 1/2 chains with Hamiltonian 
©. This problem is equivalent to a ID lattice problem 
of N sites with physical dimension (1=2x2, which can 
be easily incorporated to the framework of the numerical 
methods introduced in [TTL Il2| . Implementation of the 
scheme above is as follows. 

1) Let the state |^(0)} be the MPS with dimension 
D that better approximates the initial state |"0o) ® IV'a), 
l*(0)> = E d Sl ,..., SN =iTr(Al 1 ...A s N »)\si...s N ). Here, 
the A's are matrices whose dimension is bounded by D 
and d = 4. 2) We evolve |\&(0)) under Hamiltonian (jfJJl. 
As in [12j we take a small time step At and compute 
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|\&(Ai)) exactly, a state for which the dimension of the 
matrices A will be typically larger. Following ^1 we 
then "truncate" the matrices in an optimal way and use 
the "truncated" state to compute the next time step. 3) 
At any time t we can efficiently determine the quanti- 
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FIG. 1: Comparison of our numerical simulation of a ran- 
dom field XY spin chain with an exact calculation. We 
plot the exact results (full line) together with the numeri- 
cal results using MPS with D = 20 (circles) and D — 30 
(stars). The relative errors are plotted in the insets. Figures 
(a) and (b) show the time evolution of the averaged corre- 
lations (Ej^ + SJ^)) (blue) and «£*SfSf +1 » 
(orange). The evolution Hamiltonian is ||SJ with (a) Bo — 0, 
J/B = -2 and p(b) = 1/2 N , and (b) B = 0, J/B = -4 and 
p(b) = |ab| 2 , where a b are given by \ip a ) = 5] b Qb|b), and 
\tp a ) is the ground state of Ho with Bo /J — 1.4. For both 
figures the chain is prepared initially in the ground state of 
Ho with Bo = 0. The time step for the numerical simula- 
tions is At = O.OlJ/h. Fig. (c) shows the correlation func- 
tions ((J2 l SfSt +A + S y e S y e+A )) (blue) and <(E^fSf +A » 
(orange) averaged over all possible ground states of Hamilto- 
nian JSJ for Bo = 0, J/B — — 1 and p(h) — 1/2 as a function 
of the distance A between spins. The numerical simulation 
performs an adiabatic evolution with Hamiltonian © with 
/3(t) = Bt/T and T = lOOJ/ft. 



ties {(Oi . . . On)) as (0% ® 1 . . . On <S> 1), which can be 
efficiently computed for MPS [l^ . 

In order to test the efficiency of the simulation scheme 
above we have compared it with an exact calculation 
for the case in which Ho is an XY model Hamiltonian, 

= -JEjLi SfSf +1 + S y t S y t+1 + B Ef =1 Sf. The ex- 
act averages Q are calculated in the following way. For 
each realization of the magnetic field the evolved state (or 
ground state) of the Hamiltonian @) is computed exactly 
using fermionization techniques |2ll |22| . The quantities 
Q are then determined by averaging over the expecta- 
tion values of the 2 N states obtained. The comparison 
is shown in Fig. 1 for N = 16 and different correlations 
functions averaged over the collection of evolved states 
(Fig. 1(a) and 1(b)) and ground states (Fig. 1(c)), both 
for a uniform probability distribution (Fig. 1(a)) and a 
correlated one (Fig. 1(b)). Using MPS with D = 20 we 
obtain a very good accuracy (see insets in Fig. 1). As in 
[ITIIT^ we have two sources of error. One is the Trotter 
expansion. Here the error can be decreased by consider- 
ing smaller time steps. The other source of error is the 
truncation of the MPS to a smaller dimension. For the 
simulations we have performed this error did not grow 
much with the size of the system. Using the numerical 
scheme above we have performed simulations of the dy- 
namics of random field XY (Fig. 2) and Heisenberg (Fig. 
3) chains with N = 40. 

The scheme above can be easily extended for the nu- 




FIG. 2: Numerical simulation of the time evolution of a ran- 
dom field XY spin chain with iV = 40. We show the corre- 
lation function ({c k Ck)) as a function of time and momen- 
tum k. Here c k oc J2^sin(k£)ct, k = jf^, . . . , j^t. an d 

Ci — Il£<£' St'{Sf + iS%) are the fermionic operators given 
by the Jordan- Wigner transformation |2l|. The evolution 
Hamiltonian is @ with Ho being the XY Hamiltonian with 
Bo = 0, B/J = -4 and p(b) = 1/2 N . The initial state |^o) 
is the ground state of Ho- As the system evolves in time the 
initially sharp Fermi sea disappears. 



4 




10 20 30 40 



A 

FIG. 3: Numerical simulation of the time evolution of a ran- 
dom field Heisenberg spin chain with N = 40. We show the 
correlation function {(Yle *S!+a)} as a function of time and 
the separation A between the spins. The evolution Hamilto- 
nian is JSJ with Ho being the antiferromagnetic Heisenberg 
Hamiltonian, B/J = —4 and p(b) = l/2 N . The initial state 
\ipo) is the ground state of Hq. As the system evolves in time 
the antiferromagnetic correlations are smeared out. 



merical simulation of spin chains with random couplings. 
For example, for the case of a Heisenberg chain with 
Hamiltonian H( Ji, . . . , Jn) = J^i Ji&£ ■ Si+x, the prob- 
lem is mapped to the simulation of the time evolution of 
two spin chains under the three-body interaction Hamil- 
tonian H = H(af, . . .,$*) = J2i^e s e ■ s e+i- 

Experimental proposal. Using the ideas of the al- 
gorithm above we present an experimental scheme for 
atoms in optical lattices that could be use as a simulation 
protocol for a variety of disordered phases. We consider 
a system of atoms b (bosons or fermions) in an optical 
lattice (3D, 2D, or ID) in a certain state l^o)- We con- 
sider another system of atoms a (e.g, another spin state 
or atomic spec ies) which experiences an independent lat- 
tice potential [23. We consider a situation in which the 
lattice potentials of atoms a and b are initially shifted in 
such a way that there is no interaction between the two 
systems. We proceed as follows: 1) We prepare atoms a 
is a certain state \ip a )- This state can always be written 
in a Fock basis as \ip a ) = J2 ni ...n M &n 1 ,...n M \ n i ■ --"Af), 
where n\, . . . hm are the occupation numbers of the M 
lattice sites and the a's are certain complex coefficients. 
We then suddenly ramp up the lattice for atoms a (so 
that tunneling processes are instantaneously suppresed), 
and shift it so that the interaction with atoms b is instan- 
taneously switched on. 2) We let the composite system 
evolve. The Hamiltonian that governs the evolution is 

M 

H = H b + U ab J2nfnl (7) 
e=i 



where Hi, is the Hubbard Hamiltonian for atoms b, n£, n\ 
are the local density operators for atoms a and b, and U a b 
is the interaction coupling between atoms a and b. 3) We 
finally measure the system b. According to the algorithm 
developed above, this interacting experimental scheme is 
simulating in parallel all possible dynamics of atoms b 
under the random Hamiltonian: 

M 

H(Vx , . . . ,V M ) = H b + U ab Y, Vih\, (8) 

i=\ 

where V\ , . . . , Vm are random potential strengths that 
take values within {0,1,..., N a } with a probability dis- 
tribution given by p(V%, . . . , V M ) = \a Vl ,...v N \ 2 , N a being 
the number of a atoms. Choosing the initial state of 
atoms a in the appropriate way we can tune the prob- 
ability distribution of the random potential for atoms 
b. Changing the entanglement properties of the state in 
which atoms a are prepared we will change the correla- 
tion properties among the random local potentials. As 
well, the intensity of the disorder potential can be tuned 
by varying the interaction strength U a b (e.g., shifting the 
lattices). This allows to simulate a large variety of dis- 
ordered phases. As opposite to the classical simulation 
of randomness (e.g., with speckle lasers |24|) our quan- 
tum mechanical approach allows to simulate all possible 
evolutions of the random system in one single run of the 
experiment. Concerning measurements, note that in an 
experiment we will typically have many copies of the sys- 
tem (an array of 2D or ID systems) so that the outcomes 
will be already averaged over all copies. As an example 
of current interest let us consider the case in which atoms 
a and b are initially prepared in two independent Tonks 
states |22| with filling factor v a and v b . For this case 
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FIG. 4: Numerical simulation of the time evolution of the 
averaged density of a Tonks gas in an optical lattice with 
M — 40 sites and v = 1/2 in the random potential generated 
by another Tonks gas with v = 1/5. 
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the experimental scheme above would simulate the dy- 
namics of a Tonks gas in the presence of the random 
potential generated by another Tonks gas. Using the nu- 
merical scheme above we have simulated this situation for 
M = 40 sites and v\, = 1/2, v a = 1/5. Interestingly, the 
averaged density of atoms b shows localization of atoms b 
in the regions in which atoms a are most probably absent 
(see Fig. 4). 

In conclusion, we have presented an algorithm that al- 



lows to simulate (classical) randomness in quantum many 
body systems via a single quantum mechanical problem 
in which quantum interactions allow all possible quan- 
tum paths of the random system to occur simultaneously. 
Our scheme opens new possibilities in the numerical and 
experimental simulation of QRS. 

We thank M. A. Martm-Delgado for helpful discus- 
sions. Work supported in part by DFG, EU projects and 
Bayerische Staatsregierung. 
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